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Abstract. Semiconductor superlattices (SL) may be described by a Boltzniann- 
Poisson kinetic equation with a Bhatnagar-Gross-Krook (BGK) collision term which 
preserves charge, but not momentum or energy. Under appropriate boundary and 
voltage bias conditions, these equations exhibit time-periodic oscillations of the current 
caused by repeated nucleation and motion of charge dipole waves. Despite this clear 
nonequilibrium behavior, if we 'close' the system by attaching insulated contacts to 
the superlattice and keeping its voltage bias to zero volts, we can prove the H theorem, 
namely that a free energy of the kinetic equations is a Lyapunov functional 

((f> > 0, d^/dt < 0). Numerical simulations confirm that the free energy decays 
to its equilibrium value for a closed SL, whereas for an 'open' SL under appropriate dc 
voltage bias and contact conductivity ^(t) oscillates in time with the same frequency 
as the current self-sustained oscillations. 
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1. Introduction 

In the simplest case, semiconductor superlattices (SL) are artificial one-dimensional 
crystals comprising many periods, each consisting of two different semiconductors with 
different gaps but similar lattice constants, e.g. GaAs and AlAs. These nanostructures 
were invented by Esaki and Tsu in order to develop a device that exhibits Bloch 
oscillations [Ij. Electron transport in SL gives rise to a great variety of nonlinear 
phenomena such as pattern formation, current self-sustained oscillations and chaotic 
behavior that have been studied theoretically and experimentally [2]. These phenomena 
are best explored using balance equations for the electron density, current and energy 
densities in SL minibands and such balance equations can be derived from more 
complicated Boltzmann type equations in particular hydrodynamic limits. In practice, 
consistent derivations start from kinetic equations with BGK collision models, named 
after Bhatnagar, Gross and Krook's simplification of the Boltzmann equation for 
monoatomic gases and low density plasmas [3J. The idea behind these models is 
that the distribution function should relax to a local equilibrium distribution whose 
moments (e.g., density, average velocity and energy) are functions of the corresponding 
quantities for the instantaneous value of the distribution function. Then the Boltzmann 
collision term is replaced by a relaxation term proportional to the difference between 
the distribution function and the local equilibrium function. BGK models are popular 
because the starting point for deriving balance equations is to find a modified local 
equilibrium which is the leading order approximation in a perturbation scheme, and 
this is provided effortlessly by the BGK collision model. 

Kinetic theory descriptions of SL should explain nonequilibrium phenomena such 
as current self-sustained oscillations that are not common in the kinetic theory of gases. 
Yet the analysis of these descriptions should be relatively simple because SL are quasi 
one dimensional systems. Boltzmann type equations for SL go back to the 1970's when 
Ktitorov, Simin and Sindalovskii (KSS) |1] proposed a relaxation time kinetic equation 
that included a simple momentum-dissipating but energy-conserving collision term. 
Later Ignatov and Shashkin [5] replaced KSS's idea of relaxation to global equilibrium 
by a BGK collision term that preserved charge but dissipated energy, unlike the original 
charge, momentum and energy conserving BGK collision term [3]. More recently the 
more complete description including coupling to a Poisson equation was introduced, 
and balance equations (a nonlinear drift- diffusion equation for the electric field and the 
electron density coupled to the Poisson equation) derived by using the Chapman- Enskog 
perturbation method [6j. Among stable solutions of the balance equations there are 
self- sustained oscillations of the current through a dc voltage biased SL j6j. The kinetic 
equation has been solved numerically by a particle method providing confirmation that 
the balance equations approximate very well the field and electron density including the 
regime of self-oscillations [7]. Quantum effects [8] and electron transport in SL with two 
populated minibands have also been considered [HI [10] . 

In this work, we find the Lyapunov functional for a BGK-Poisson kinetic equation 
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describing a SL for idealized boundary conditions (infinite SL, periodic boundary 
conditions or finite SL with insulated contacts at zero voltage bias). This functional 
has the form of a free energy and the corresponding entropy production contains 
a generalized force which is proportional to the difference between distribution and 
local equilibrium distribution functions. It is remarkable that formulas of irreversible 
thermodynamics hold in a regime that is very far from equilibrium. We approximate the 
free energy functional using the leading order approximation of the Chapman-Enskog 
method and calculate it numerically in two cases: insulated contacts at zero voltage bias 
{closed system) and contacts having finite conductivity at nonzero voltage bias {open 
system) in a regime where there are self-sustained oscillations of the current through the 
SL. In the first case, we have checked that the free energy indeed decreases monotonically 
towards its equilibrium value, whereas for the case of the open system (nonzero voltage 
bias and nonzero contact conductivity) the free energy oscillates in time following the 
evolution of the current and the electric field inside the SL. 

The rest of this paper is as follows. Section [2] presents the model kinetic equation 
and a generalized drift-diffusion system that can be derived from it by means of the 
Chapman-Enskog method [6J. In Section [3l we derive an entropy density functional 
from a given form of the local equilibrium distribution, and then find a free energy 
density from which we prove the H theorem provided the contacts at the SL ends are 
insulating and there is a zero voltage bias between them. In such conditions, any initial 
condition evolves towards the globally stable equilibrium and the free energy decreases 
monotonically until it reaches its equilibrium value. Of course under realistic boundary 
conditions of charge injecting and collecting contacts and nonzero voltage bias, the SL is 
far from equilibrium and, in fact, self-sustained oscillations of the current due to periodic 
recycling and motion of charge dipole domains are possible stable solutions [2], [7]. In 
this case, the free energy is no longer a Lyapunov functional and it oscillates periodically 
in time in the regime of current self-oscillations. In Section Hj these different regimes 
and behaviors of the free energy are confirmed by solving numerically the generalized 
drift-diffusion system, reconstructing the approximate distribution function provided by 
the Chapman-Enskog method and finding the free energy as a function of time. The 
last section contains our conclusions. 

2. BPBGK kinetic equation and generalized drift-diffusion balance equation 

In this Section, we present the Boltzmann-Poisson kinetic description of a one-miniband 
SL and recall the drift- diffusion balance equation that can be derived therefrom by means 
of the Chapman-Enskog method [6]. 
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2.1. Kinetic equation and local equilibrium distribution 



The Boltzmann-Poisson-Bhatnagar-Gross-Krook (BPBGK) system for ID electron 
transport in the lowest miniband of a strongly coupled SL is [Sj: 



df_ ld£df_ e dW df _ 
dt h dk dx h dx dk 
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Here /, Af = [/(x, k, t-f{x, -k, t)]/2, f^^{k; n), n, Nd, S{k), ds, dw. I = dn+dw, W, 
e, m*, kB, T, T, z/e„, Uimp and — e < are the one-particle distribution function, its odd 
part (in k), the ID local equilibrium distribution function, the 2D electron density, the 
2D doping density, the miniband dispersion relation {S{k) is even: £{—k) = S{k)), the 
barrier width, the well width, the SL period, the electric potential, the SL permittivity, 
the effective mass of the electron in the lateral directions, the Boltzmann constant, 
the lattice temperature, the energy broadening of the equilibrium distribution due to 
collisions [11] (page 28), the constant frequency of the inelastic collisions responsible 
for energy relaxation, the constant frequency of the elastic impurity collisions and the 
electron charge, respectively. The chemical potential is a function of the electron 
density n that can be obtained by inserting (jl]) into ([3]) and solving for /i = fj,{n). Thus 
/i is a functional of / which we may write as 

Integrating ([1]) over k and using given by ([3]), we find the charge continuity 
equation: 

dJn 



e dn 
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0, 



where J„(x,t) is the electron current density: 

f7r/« 
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1 dS 



v{k) f{x,k,t)dk, v{k) = - — 



h dk 
we obtain 
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and v{k) is the electron group velocity. From ([2]) and ([5]) 

OF dW 
e^ + J„^m, F^-, (7) 

where —F is the electric field, J{t) is the total current density and ([7]) is a form of 
Ampere's law. The main idea behind the BGK collision model is to substitute the 
linear, quadratic or quartic Boltzmann collision terms, which are nonlocal in k, by a 
nonlinear collision term, Q[f], that is local and preserves charge continuity, as in ([5]). 
It is convenient to use the following tight binding dispersion relation for the miniband 
with lowest energy: 
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Eq. dH]) contains the first two harmonics in a Fourier series of the periodic function £{k). 
dH]) approximates the band dispersion relation of any crystal provided the band width 
A is small compared to the energy difference between bands (band gaps) [12]. Besides 
being a reasonable approximation for a SL with well-separated minibands, ([8]) produces 
simple analytic expressions for the balance equations and is often used in comparisons 
between theory and experiments [2]. 

The system ([1])-(11]) is the semiclassical limit of the Wigner system of equations 
considered in Section 2 of Ref. [9]. Note that the ID Fermi-Dirac local equilibrium 
distribution (HI) can be written as 



f^{k,n)=^{E[f]\ E[f] = ^^^, (9) 

7(i?) is a decreasing function that takes on positive values for real values of E. If 
broadening due to scattering is negligible, F — )■ 0, and ffTOj) becomes 

^(E) = -^Hl + e-^). (11) 

In the Boltzmann limit, E oo, and we have 

, , m*kBT p , , 

= (12' 

This local equilibrium distribution was used by Ignatov and Shashkin [5] to analyze 
particular solutions of the BPBGK system for constant electric field. 

2.2. Chapman- Enskog perturbation method and drift- diffusion system 

In the hyperbolic limit where collisions and field-dependent terms in the BPBGK system 
dominate all other terms, it is possible to derive a generalized drift- diffusion equation 
for the electric field, —F = —dW/dx. Collision and field-dependent terms in ([1]) are 
of the same order provided e[F]l/{h) = Uen, where [F] gives the order of magnitude of 
the field. Then [F] = hvenKd). The Poisson equation implies that the distance [x] 
over which field varies an amount [F] is proportional to s[F] divided by cNd/I, so that 
[x] = e[F]l / (cNd) = ehuen/ (c.'^Nd). Let vm be the order of magnitude of the electron 
velocity (see below). The reciprocal of the electron residence time is vm/[x] and the 
condition that the field-dependent and collision terms dominate all others in ([T]) is: 

Provided this condition and ([8]) hold, the terms df/dt and v{k) df /dx, with hv{k) = 
S'{k) = dS/dk, are both of order 5^1 compared to ehr^dW / dx df / dk and Q[f] in 
([T]). Ignoring these small terms, ^ can be approximated by 

£/W(A;; F) ^ {^-^^ + u,^ + u,^,A^ f<'\k- F) = u^J^'' . (14) 



< ^ 5 = — < 1. (13) 
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The solution of this equation, 

oo 
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is the Fermi-Dirac local equilibrium modified by a high electric field. The field-dependent 
local equilibrium /"^^^ is the leading order term of a Chapman-Enskog expansion 

/~/W(A;;F) + /«(A;;F) + /(2)(fc;F), (17) 

e^ + J^''\F]+J^^\F]^J{t). (18) 

(JTTj) becomes a power series in the small dimensionless parameter 5 of ( |T3i) if we 
nondimensionalize the BPBGK system using the scales J-", x = a;/[x], A; = kl, 
t = VMt/lx], etc. The corrections f^"^\ m = 1,2, solve the equations 



(1) 



5/(0) 



(2) 



9t 



dt 



-v{k) 



- t;(A;) 



5/(o) 
dx 



d 



f 



(0) 



(19) 



(20) 



dx dt' 

in which the subscripts and 1 in their right hand side (RHS) mean that edF/dt is 
replaced by J — J'^^^[F] and by —J''^^^[F], respectively. The conditions that /(^) and 
Z*-^^ be bounded and 2ti/1 periodic in k determine the functionals J'^^^[F] [9J. When 
these functionals are inserted in ( TT8l) . we obtain the following generalized drift-diffusion 
system for the field P, '9] 
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In these equations, fi = M provided n = N^. Equation ( 12T|) is a form of Ampere's law 
establishing that the sum of the electron current density and Maxwell's displacement 
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current equals the total current density. Numerical solution of both the kinetic system 
([T])-(j4]) and of the drift-diffusion system fl2T]) - fl28|) show that the later is a very good 
approximation to the former for realistic parameter ranges [7j. 



3. H theorem and entropy for ideal boundary conditions 

In this section, we find a Lyapunov function for the BPBGK kinetic system assuming 
that the SL is infinite or that it is finite but it has ideal boundary conditions at the 
contact regions (zero voltage bias and either periodic boundary conditions or insulating 
contacts with zero conductivity). The Lyapunov functional has the form of a free energy 
and, for a local equilibrium of Boltzmann type, it produces formulas that are very close 
to those of equilibrium thermodynamics. 



3.1. Derivation of a free energy density 

As in the case of the usual BGK equation in the kinetic theory of gases, let us start 
by finding an entropy functional [13]. Suppose we have a local equilibrium 7(-E), with 
E = E[f], which is a non-negative decreasing function for real values of E. 'y{E) can 
be found by maximizing the entropy functional 

S{x,t) = — dk ds-f-\s) (29) 

2vr J^n/i Jo 

(where 7~^(s) is the inverse function of 7(-E); S has units of 3D density), subject to the 
conditions 

dkf{x,k,t)=n{x,t), (30) 

n/l 

I r'^ 

dkS{k) f{x,k,t) = S{x,t). (31) 

n/l 
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In fact, taking the extremal of 

S(x,t) + i^^ r'\kf{xXt)-^^ r" dk8{k)f{xXt). (32) 

we obtain 'j''^{f^^) = f3S{k) — fi. Then 

r^ = 7(/3^(fc)-/i). (33) 

The Lagrange multipliers (3 and fi have to be calculated by inserting ( 133|) in (!30l) - (!3T!) 
and solving these equations for f3{x,t) and fi{x,t) in terms of n{x,t) and S{x,t). The 
entropy density S was introduced by Aoki et al [14J for a BGK kinetic equation with 
charge and energy conserving collision terms representing infinitely extended materials 
with parabolic bands. If we set (3 = 1/(A;bT), fi = fi/^ksT), and, by an appropriate 
choice of fi{x,t) with constant T, we impose that ( j30l) be still satisfied but not ( l3T]) . 
f^D gjygj^ ]-,y ( |33|) becomes the same function of (|9l)- (|T0l) . 
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To find the Lyapunov functional, we define a free energy density 

r]{x,t) = v{x,t) - S{x,t)l = v{x,t) / dk rfs7"^(s), (34) 

27r J^^/i Jo 

wliere v{x,t) is a functional of the distribution function to be specified later. Using 
periodicity in k, fl29|l and flM|) yield 
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The first term in the RHS of this expression is nonpositive because the integrand is 
proportional to the derivative of the decreasing function 7"^ at some intermediate 
point (mean value theorem). The same argument shows that the third term in 
the RHS is also nonpositive. This term has been written as indicated using that 
j:i)^g{k)Afdk = j:i)^AgAfdk for any function g{k). The fourth term is the x 
derivative of some function. The second term in the RHS of (|35l) is zero for energy 
conserving collision terms [13] but this is not the case for the BPBGK system. Let us 
rewrite this term by using 
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after integration by parts. The last term in the RHS of (1361) is 
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Then §^ and §3) yield 
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We have omitted the dependance on x, t in the second term of the RHS of (l39l) to 
simphfy this formula. Selecting now 
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The free energy production cr^ in (jUj) is zero only if / = 7(i?[/]) = f^^ (which is even 
in k). Moreover 7~^(s) is decreasing from infinity at s = and therefore it vanishes at 
a unique value s = Sq > 0. This value yields the maximum of the integral 'y~^{s)ds. 
Let us select 

C= rdsr\s), (45) 



then the free energy density becomes 
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3.2. Boltzmann limit and physical meaning of ri[x,t) 



In the Boltzmann limit, ( 12^ yields 7 ^(s) = — \n{s/Q), with Q = m*kBT / [iifi?) 
- Jol'Hs) ds = f ln(//g) ~ f, so = C = Q because of (05]), and (gSD gives 
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The first two terms in the RHS of ( 1471) are the internal energy density (including the 
contribution from the electric field) and, except for an additive constant, the two last 
terms can be written as —Ts{x,t), the temperature times the entropy density if we 
define 

Thus ksTr^/l is the free energy density. It is remarkable that the classical expression 
for the free energy holds exactly for out of equilibrium states described by the BPBGK 
equation in the Boltzmann limit. For other local equilibrium distributions, the entropy 
functional ( |29l) does not yield the recognizable extensive form of entropy. In the 
Boltzmann limit, ([3]) gives 
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which is the usual logarithmic dependence of the chemical potential on the electron 
density. The explicit formula in f l5U]) . where Io{x) is a modified Bessel function, is found 
for the tight binding dispersion relation (|8]). If the distribution function is the local 
equihbrium ( !T2|) . the free energy density (H7|) becomes 
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3.3. Boundary conditions and free energy as a Lyapunov functional 

For an infinitely long SL under zero voltage bias conditions or for a finite SL with 
periodic boundary conditions, the free energy flux (H3l) obeys J^(L,t) = J^(0,t), and 
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and ( HBj) imply that the total free energy, 



, kBTA , , , 



(52) 



is a Lyapunov functional of the BPBGK system ($(t) > 0, < 0; A is the area 

of the SL cross section). Then the local equilibrium with n = Njy and F = is the 
globally stable stationary solution. These idealized boundary conditions do not hold for 
a finite SL under voltage bias and with realistic boundary conditions [7]. We have 



which no longer has a definite sign. 

What are the boundary conditions used to solve the BPBGK system? For a finite 
SL of length L = (N + l)l {N is the number of SL periods) and a tight binding dispersion 
relation, appropriate boundary conditions are [7] 



These relations are of the forms: J„(0,t) = jo(-^(0,t)) and J„(L,t) = 
jL{F{L,t))n{L,t)/ND in which the contact currents jo{F) and jiiF) have been 
linearized according to Ohm's law (with contact conductivities (Tq and ai, respectively). 
A derivation of these relations from quantum theory can be found in [15] . They are often 
used in the literature; cf the review [2]. We prove below that J^(L,t) — J^(0,t) > 
provided ctq = o"l = (insulated contacts) and V = Q (zero voltage bias). Then the total 
free energy, $(t), is a Lyapunov functional of the BPBGK system and the distribution 
function tends to equilibrium which is the unique globally stable stationary solution of 
the system. 

Proof that Jr^{L, t) — Jrf{0, t) > for zero contact conductivity and zero voltage bias. 
(Adapted from Chapter 8 of [IS]). In (l44l) . the function 



is convex for non- negative / and it reaches a minimum at / = / . Then Jensen's 
inequality implies that 




(53) 





(57) 




(58) 
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at X = L, where 

m k) = = vik')mk) 

H/M'^)\f^'Kn)dK Jl^,v{K)f(^)iK)dn 

with A;' > and A; < 0. We have omitted the dependence on x = L and t in the 
distribution functions. In the left hand side of 

fTr// 

n{k',k)f{k')dk' = f{k), (60) 



according to with = 0. We now substitute (pi]) in ( 155]) . multiply the result by 
= —v{k), integrate over k in — vr/Z < A; < 0, and use that — J_^^^v{k)Q{k' , k)dk = 
v{k') due to dsn])- The result is 

v{k)C{f{k))dk< / v{k)C{f{k))dk, therefore 

n/l Jo 

n/l 

v{k)C{f{k))dk>0. (61) 

-tt/I 

Similarly, we can prove that j^'^jiv{k)C{f{k))dk < at x = if ctq = 0. Since 
= for insulating contacts, (H3|) implies Jr^{L,t) — J^(0,i(:) > —JlV/{kBT) = for 
V = (zero voltage bias). Then and (gl]) yield d^/dt < 0. 



4. Numerical results 



In order to calculate numerically the Lyapunov functional $(t), we have first obtained 
F{x,t), n{x,t) and Jit) by solving f l2T|) - fl28|) under appropriate dc voltage bias and the 
boundary conditions ( 156|) . From equation ( IT5|) we have calculated the leading order 
approximation of / ~ /*^°^ using the Boltzmann local equilibrium distribution ( IT2|) (in 
this case = ^ I j{ 2KBT ^ / ^^^ 2KbT ^ ^^^^^ -^i(^) is a modified Bessel function). Then 
we have obtained the free energy density rj from fHTl) . Finally, we have calculated the 
Lyapunov functional ^(t) from equation fl52|) . 

In our numerical simulations, we have used the following typical parameters |2]|: 
(ifi = 1-5 nm, dw = 9 nm, I = ds + dw = 10.5 nm, A^^d = 2.5 x 10^° cm~^, 
i/en = 9 X 10^2 Hz, Uimp = 0, = 80, T = 5 K, A = 2.6 meV, £ = 12.85^0, 
A = 100 yum^. We have selected the following units to present our results graphically: 
Fm = huen/{el) = 5.64 kV/cm, xq = eFMl/ieNo) = 16.83 nm, vm = M/{2h) = 20.46 
km/s, Jo = eNDVM/l = 7.8 kA/cm^, to = Xo/vm = 0.82 ps, r/o = ^ffgj = 2.76 x lO^^ 
cm~^, $0 = kBTArjoXo/l = 190.27 eV. The chosen unit rjQ reflects the fact that the 
most important contribution to the free energy density is the electrostatic one. 



We have studied the following cases: 

Zero contact conductivity and zero voltage bias: We see in figure [1] that the 
Lyapunov functional $(t) tends asymptotically to the local equilibrium ^^{t) = 
^^f^ lo V^{x, t) dx, as expected. 
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Figure 1. Time evolution of the Lyapunov functional ^{t) for zero contact 
conductivity and zero voltage bias. Here $0 = 190.27 eV, to — 0.82 ps. 
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Figure 2. (a) Total current density vs time, (b) Electric field profiles at several times 
during one oscillation period. Here Jq = 7.8 kA/cm^, to — 0.82 ps, Fm — 5.64 kV/cm, 
xq — 16.83 nm. 
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Figure 3. (a) Time evolution of the Lyapunov functional during self sustained 
current oscillations in a SL with nonzero contact conductivity biased at 1.7 V. (b) 
Same for $(t) - $^(t). 
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• Realistic contact conditions: For an applied voltage bias V = 1.7 Y and contact 
conductivities ctq = 1.38fi~^m~^ and (Jl = 0.69 f2~^m~^, the total current exhibits 
self sustained oscillations due to the repeated nucleation of dipole waves at the 
injecting contact and their motion towards the collector [2], as shown in Figure |2j 
Figure [3] shows that free energy oscillates at the same frequency and, obviously, 
it is no longer a Lyapunov functional. Global equilibrium is precluded by injection 
and depletion of carriers at the contacts and the nonzero voltage bias. Figure [3]^b) 
shows that the free energy differs from that of local equilibrium whereas a compar- 
ison to Figure [3](a) indicates that has the same order of magnitude as $. 

5. Conclusions 

Electron transport in strongly coupled miniband superlattices is described by a BPBGK 
kinetic system of equations. For this system, we have found a free energy density 
ri{x,t) and its corresponding free energy We have proved that the free energy is 

a Lyapunov functional ($ > 0, d^/dt < 0) for idealized boundary conditions in the SL 
(zero voltage bias and either zero contact conductivity or periodic boundary conditions). 
Thus for a closed SL with insulating contacts and zero voltage bias, the distribution 
function tends to equilibrium which is the unique globally stable stationary solution of 
the BPBGK system. In the Boltzmann limit, it is remarkable that the expression of the 
free energy density rj{x, t) is similar to that of equilibrium thermodynamics, including the 
internal energy density minus a T s{x,t) term (temperature times the entropy density). 
However, an open SL (nonzero voltage and nonzero conductivity at contacts) may be 
very far from equilibrium and display self-sustained current oscillations due to repeated 
nucleation and motion of charge dipole waves as in Figure El In this case, c/$ / dt has 
no definite sign. Numerical solutions of the BPBGK system of equations illustrate the 
monotonic decay of the free energy of a closed SL towards its equilibrium value and 
show that, during self-oscillations in an open SL, $(t) oscillates stably and periodically 
in time at the same frequency as the total current. 

Our results are based on selecting an entropy density based on the local equilibrium 
distribution 7 (£'[/]) by using the maximum entropy principle and therefore they can be 
applied to other forms of 7(-E) (decreasing with E, compact support [H]) or to higher 
dimensional models of similar type (collision operator preserving charge but not energy 
or momentum). 
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